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ABSTRACT 
The  semi-implicit  method  for  the  MHD  equations   is  extended  to 
grid-point  models  and  for   implementation  with  a  more  accurate  time 
advance  scheme. 
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I.   INTRODUCTION 

In  many  nonlinear  magnetohydrodynamics  (MHD)  problems  it  is 
important  to  retain  most  of  the  physics  contained  in  the  resistive  MHD 
equations.  Approximations,  such  as  the  assumption  of  incompressibility 
or  the  use  of  an  inverse  aspect  ratio  expansion,  can  sometimes  lead  to 
erroneous  results.  A  major  difficulty  in  solving  the  full  MHD 
equations  has  been  the  fact  that  when  features  like  compressibility  are 
retained,  parasitic  waves  (e.g.,  short  wavelength  magnetosonic  modes) 
will  also  be  present.  These  waves  generally  have  little  effect  on  the 
long-time  solution,  however,  they  can  force  explicit  codes  to  use  time 
steps  governed  by  such  a  restrictive  Courant-Fr iedrichs-Lewy  (CFL) 
condition  that  economical  long  time  scale  computations  are  impossible. 
Semi-implicit  methods  have  been  used  successfully  to  solve  such 
problems  in  other  fields,  particularly  meteorology.  In  plasma  physics, 
the  semi-implicit  MHD  method  '-''  provides  a  simple  means  of 
eliminating  CFL  constraints  due  to  parasitic  waves,  while  still 
treating  the  long  time  evolution  of  the  resistive  MHD  equations 
properly.  This  paper  describes  two  straightforward  enhancements  of  the 
basic  semi-implicit  MHD  method.  Section  II  describes  an  implementation 
of  the  method  in  which  the  semi-implicit  terms  are  treated  with  fourth 
order  accuracy  in  time.  Section  III  describes  how  to  implement  the 
method  in  a  grid-point  model. 
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II.   TIME  ADVANCE 

In  the  standard  semi-implicit  MHD  methods  '2'  the  semi-implicit 
terms  enter  with  second  order  accuracy  in  time.  It  is  possible  to 
reduce  further  the  error  arising  from  the  semi-implicit  terms  by  making 
an  additional  pass  on  the  momentum  equation  so  that  the  semi-implicit 
terms  are  treated  effectively  at  the  same  time  level.  The  error  due  to 
these  terms  is  then  reduced  to  fourth  order  in  time. 

For  simplicity,  consider  the  MHD  equations  with  zero  pressure, 
uniform  density  (p  =  1.0),  and  neglecting  advection.  The  MHD  equations 
then  become 


3V      ■*    ^ 

_  =  (VxB)  =<  B  (1) 


4r  =  ^  ^  (VxB)  .  (2) 

ot 


A  simple  semi-implicit  advance  for  these  equations  is  the  leap-frog 
4 


scheme i 


1       1 
B   2  =  B   2  +  AtV  X  (v"xB")  (3) 


1        1 
P  P^        ^n+ 

oi  ol 


V"""^  -  A2,(At)2v2v"^^  =  v"  -  A^  (At)2v2v"  +  AtCVxs"^  2)  x  b"^^  2   (4) 


The  simplest  form,  the  Laplacian,  has  been  chosen  for  the  form  of  the 

p 

semi-implicit  operator.    The  coefficient,   A   ,   is  the  sum  of  the 

squares  of  the  components  of  A^  perpendicular  to  V.  As  an  example,  in 

?     ?     ?  '* 

the  V   advance  A  ,  =  A^  +  k^    .      It  is  important  for  accuracy  that  A^ 

■^  Q  '-'''  '-'J  O 
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be  treated  as  a  vector  rather  than  a  scalar   in  cases   in  which  one 

component  of  the  magnetic  field  dominates,  such  as  the  toroidal  field 

-> 
in  a  tokamak.   The  magnitudes  of  the  components  of  Aq  are  chosen  from 

stability  considerations. 

For   linear   waves  the  algorithm  of  Eqs.    (3)   and  (H)      is 

p 
unconditionally  stable  provided  A    is  sufficiently  large.    In  one 

°^ 

?    1   ? 
dimension  this  requires  simply  that  A   >  -  B^.   The  algorithm  is 

ol   ^  1 
non-dissipative  and  second  order  accurate   in  time.    As  with  any 

leap-frog  method,   dissipation  may  be  added  if  it  is  desired.   The 

Alfven  waves  produced  by  Eqs.   (3)  and  (4)  have  the  dispersion  relation 

sm'^  [^J  = .  (5) 

2     1  +  A2(At)2 


The  solution  may  be  made  more  accurate  by  adding  one  more  pass  to  the 
momentum  equation  solution,  giving  the  scheme 


1       1 
-»n+  —   ->n-  —         ->■->• 

B   2  =  B   2  +  AtV  X  (v"xB")  (6) 


* 


1        1 
n+  -    ->-n+  — 


A^  (At)^V^V   =  v"  -  A"^  (At)^V^v"  +  At(VxB    2)  x  B  2      (7) 
ol                 ol 

^         ->•  +  ^  ^  +  ^ 

V"*^  -  A^  (At)2v2v"^^  =  v"  -  A^  (At)2v2v*  +  AtCVxe"   2)  x  b"   2   (8) 


o 


1  ol 


The  error  due  to  the  semi-implicit  terms  is  fourth  order  in  time  in 
this  algorithm.  The  dispersion  relation  becomes 


1  +  2A^  (At)^ 

sm^  l-r-J  =  k-^V^  _ ± (9) 

2  1  +  2a2  (At)2  +  A\(At)^ 

4         4 

2     1   2 
Unconditional  stability  now  requires  A   >  _  B  .   The  two  dispersion 

4   2  [ 
relations  are  plotted  in  Fig.   1  using  their  respective  minimum  values 

2 
of  A   .   Even  though  both  schemes  are  second  order  accurate,  the  fourth 

4 

order  treatment  of  the  semi-implicit  terms  provides  a  substantial 
improvement.  The  proper  Alfven  wave  frequency  is  recovered  within  one 
percent  even  for   time   steps   larger  than   w  At  =  O.H.      The  second 

Pi. 

algorithm  does  require  more  computation  time  per  time  steps  than  the 
standard  one.  For  the  full  MHD  equations  in  a  three-dimensional 
spectral  code,  however,  the  increase  is  normally  less  than  20^.  The 
small  increase  in  time  is  due  to  the  fact  that  the  computation  time  for 
a  multi-mode  case  is  dominated  by  the  convolutions  (for  either  fast 
Fourier  transforms  or  direct  summation).  The  additional  equation  [Eq. 
(8)]  only  adds  one  additional  cross  product  to  the  algorithm,  which  may 
alternatively  be  stored  from  the  prediction  to  give  a  negligible  time 
increase.  The  increase  in  storage  requirements  will  also  be  less  than 
205^.  The  second  algorithm  is  preferable  because  its  higher  accuracy 
can  be  obtained  at  a  relatively  small  additional  cost.  In  Ref.  3  the 
same  procedure  as  in  Eqs.  (6)-(8)  was  used  in  a  predictor-corrector 
algorithm  to  obtain  accurate  results  at  large  time  steps. 
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FIGURE  1 


Dispersion  relation  for  the  semi-implicit  method  with 
second  (solid)  and  fourth  (dashed)  order  semi-implicit 
terms. 
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III.   GRID-POINT  MODELS 

The  most  obvious  application  of  the  semi-implicit  method  is  in  a 
spectral  model.  In  these  cases  the  computation  time  is  dominated  by 
performing  convolutions  and  the  additional  expense  of  the  semi-implicit 
algorithm  consists  only  of  solving  some  relatively  inexpensive  block 
tridiagonal  matrix  systems  with  at  most  S^S  blocks.  The  semi-implicit 
method,  however,  is  also  well  suited  for  implementation  in  grid-point 
models.  One  way  of  solving  the  semi-implicit  equation  is  with  a  direct 
solver,  provided  sufficient  storage  is  available.  Since  the 
semi-implicit  coefficients  are  only  approximations  to  the  magnetic 
field  values,  the  LU  decomposition  does  not  have  to  be  regularly 
updated.  It  may  be  stored  once  and  then  retained  as  long  as  the 
coefficients  satisfy  the  stability  conditions.  A  simpler  way  to  use 
the  semi-implicit  method  in  a  grid-point  model   is  with  a  factorized 

operator.    Factorized   operators   have   been   used   recently  in 

5 
semi-implicit  meteorological  calculations  by  Tanguey  and  Robert.   For 

the  MHD  equations   this  may  be  done  very  simply  when  the  Laplacian  is 

used  as  the  semi-implicit  operator. 

Consider   Eq.      (4)    in  two-dimensional   Cartesian   coordinates. 


[1    -a2    (At)2(4.4jK-1 

=    [l    -    a2    (At)2[J_   +   A_)]V"    +    At(VxB"^    2)    X    b"      2    .       (10) 
ol  3x2        9y2 


The  semi-implicit  operator  may  be  approximated  by 


1  -  a2  (At)2[4  ^  4^  ^  A^At)^  _|1_  =  1  -  AVAt)2[4  .  4)   (11) 


as  in  alternating  direction  implicit  (ADl)  methods.  In  order  to  use 
this  form  in  the  serai-implicit  equations  one  can  use  a  factorization 
like  that  used  by  Beam  and  Warming,   so  that  Eq.   (10)  becomes 

V*  =  [1  -  a2  (At)2  X^]    [1  -  a2  (At)2  iljv"  (12) 

(1  -  A^  (At)2  ^]V    =  V   +  At(VxB    2)  X  B    2  (13) 


(1    -   a2    (At)2ii]>l    =    ?' 
4  3y2 


(U) 


The  stability  limits  for  this  algorithm  are  identical  to  those  for  the 
unf actor ized  operator.  An  additional  source  of  error  has  been 
introduced  due  to  the  factorization.  For  reasonable  time  steps, 
however,  this  factorization  error  is  negligible.  In  Fig,  2  the 
dispersion  relation  for  Alfven  waves  propagating  at  45°  to  the 
coordinate  axes  (the  worst  case  for  factorization  error)  is  pJjtted 
using  both  of  the  leap-frog  algorithms  of  Sec,  II.  For  w^At  <  O.^J  the 
results  are  virtually  identical  to  those  using  the  unfactorized 
operator.  The  algorithm  is  inexpensive  as  it  requires  only  two  simple 
tridiagonal  solutions  in  two  dimensions  (three  in  three  dimensions). 
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A 


Figure  2 


Dispersion  relation  for  the  factorized  method  with 
second  (solid)  and  fourth  (dashed)  order  semi-implicit 
terms. 
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IV .   SUMMARY 

The  semi-implicit  MHD  method  provides  an  efficient  means  of 
solving  long  time  scale  MHD  problems.  It  is  straightforward  to 
implement  the  method  in  a  grid-point  model  as  well  as  in  a  spectral 
model.  Accuracy  may  be  improved  by  making  the  semi-implicit  terms 
fourth  order  accurate  in  time  with  a  relatively  small  additional 
expense  in  terms  of  computation  time. 
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